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Abstract. Based on the scalar Helmholtz equation and the finite-difference 
approximation, we formulate a matrix eigenvalue problem for the calculation of 
propagation constants, P(u)), in micro-structured optical fibres. The method is 
applied to index-guiding fibres as well as air-core photonic bandgap fibres, and in 
both cases qualitatively correct results are found. The strength of this approach 
lies in its very simple numerical implementation and its ability to find eigenmodes 
at a specific eigenvalue, which is of great interest, when modelling defect modes 
in photonic bandgap fibres. 
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1. Introduction 

Since the first resuhs on photonic crystal fibres (PCF) many exciting phenomena 
have been reported (for a recent review see e.g. Ref. and references therein). From 
the very start, a great emphasis has been on the modelling of the optical properties 
and frameworks of high complexity have been developed. In this work, we develop a 
"poor man's approach" which allows for easy implementation and calculation of the 
propagation constant, /3(w), in arbitrarily microstructured fibres. Most approaches 
start from the fundamental vectorial wave equation for an isotropic dielectric medium 

V X -i^V X H(r) = k^H(r) (1) 
e(r) 

where k = u/c and e(r) is the dielectric function, which may be frequency dependent. 
For a fibre geometry with e(r) — e{x,y), i.e., translational invariance along the z- 
direction, we look for solutions of the form H{r) = h{x,y)e^'^^ . Substituting this 
ansatz into the wave equation results in an eigenvalue problem, which determines 
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Figure 1. Contour plot of Eq. for u)^ = ^ [/3^ + -fsT^] shown by the dashed 
hne. For the region enclosed by the solid line, the relative error of the finite- 
difference approximation is generally less than 9%, and the relative error is zero 
at the origo {kx = Ky = 0) where u)'^ = ^0"^ . 



the dispersion ijj{(3). In the hterature, it is often emphasized that in general a 
ftilly-vectorial approach is required for quantitative modehing of micro-structured 
fibres. Several fuUy-vectorial approaches have been reported including plane-wave 
methods 13 multi-pole methods (SJ E], localised- function methods |7| and 
finite-element methods The complicated implementation is common to all these 
methods. In this work, we present a "poor man's approach" for calculating the 
propagation constant, in arbitrary dielectric structures. Despite its obvious 

shortcomings, we believe it is far more useful for more qualitative studies as well as 
for teaching of the physics of micro-structured optical fibres at the introductory level. 



2. Theory 

2.1. The scalar Helmholtz wave-equation 

We start from the scalar Helmholtz wave-equation 

{£2+-^+ £(2;, y)fc') y) = P'^ix, y) (2) 

with ^ being the scalar field. In this approximation, we have neglected a logarithmic 
derivative of the dielectric function at the interfaces between e.g. Si02 and air. We also 
note that polarisation effects and/or degeneracies may only be fully revealed by a fuUy- 
vectorial approach. The strength of the scalar approach is that it is straight forward to 
implement and thus serves as an excellent starting point for researchers and students 
entering the field of micro-structured optical fibres. Furthermore, the Helmholtz 
equation ^ allows for an easy incorporation of effects like material dispersion, e(w). 
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2.2. The finite- difference approximation 

Equation ^ constitutes an eigenvalue problem from which /3(ijj) can be calculated. 
We take the probably most simple approach based on a finite difference approximation 
of the Laplacian. For a quadratic grid with j labeling the grid points with spacing a 
we e.g. get the symmetrized representation of the Laplacian 

^f[x = ja] « ^ (/[(j + + f[{j - l)a] - 2f[ja]) (3) 

corresponding to nearest-neighbour coupling of the grid points. We can thus rewrite 
Eq. as a matrix eigenvalue problem 

0* = (3^<b (4) 

with 

j = i 

j, i nearest neighbours (5) 
otherwise 

where K — 1/a. For the dielectric fmiction we have ej — e{xj,yj) where {xj,yj) is 
the coordinates of the jth lattice point. The numerical task thus consists of finding 
eigenvalues of the matrix 0, which is easily done using standard numerical libraries. 
The matrix is highly sparse, symmetric and when the dielectric function is real, it is 
even Hermitian. The numerical accuracy of the approximation is of course increased 
by decreasing the lattice constant. 

2.3. The homogeneous case 

In order to estimate the required size of a, we first consider the homogeneous case 
where ej = e. This problem is well-known from solid state theory; it can be 
diagonalized by a plane wave ansatz, which results in the usual cosine-band result 




^ — [P'^ + 2K'^ {2 ~ COS Kr^a - COS Kyo)] . (6) 

This result also has the correct asymptotic behaviour of the homogeneous-space 
solution 



and by choosing a sufficiently small, the numerical discretisation is a good 
approximation to the solution of the exact problem. Because of the discretisation 
procedure, Eq. 10 has a finite band- width of ^ ^ 8K^, i.e., 

max{2K^ {2 — cos K^a ~ cos Kyo)} (8) 
— min{2i4r^(2 — cos KxO ~ cos Hyo)} = 8K^. 
This means that only frequencies satisfying 

<CJ^<-\P^ + 8K^] (9) 
e e 

can be accounted for by the discretisation procedure. In the low-frequency regime, 
the relative error of the finite-difference approximation becomes small, and typically 
we should be limiting ourselves to e.g., 

Lo'<j [P^ + K^]^a< {eu?l^ - fi^Y^'^ (10) 
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where the relative error of the finite-difference approximation is less than 9%, see 
Fig. n For higher frequencies, the finite-difference procedure becomes an inaccurate 
approximation to the exact problem, because of artificial band-structure effects. 

2.4. The boundary problem 

In principle is infinite and for the implementation we obviously need to truncate 
the matrix. This truncation may affect the accuracy of the calculation. However, we 
are also faced with the problem of deciding how the finite-difference representation of 
the differential operators (in our case the Laplacian) are represented on the boundary 
of the calculation domain. This problem arises because calculation of finite differences 
on the boundary requires the use of points that fall outside the calculation domain. 
Thus, we have to determine a proper way of how these non-existing points should be 
treated. 

In the definition of the ©-matrix, we have simply chosen to neglect the points that 
fall outside the calculation domain. This is done in order to limit the complexity of the 
0-matrix, and thereby keep the formulation of the problem as simple as possible. The 
consequence of this simplification is that we impose the condition that the field has to 
be zero on the boundaries of the calculation domain. Naturally, this assumption will 
have an effect on the accuracy of the calculation, but the better the field is confined 
within the calculation domain, the better the truncated problem resembles the correct 
solution, since a zero amplitude on the boundary becomes a reasonable approximation 
in this case. 

3. Modelling of Photonic Crystal Fibres 

3.1. Numerical Implementation 

Once the theory of the finite difference approximation has been established, the task 
of finding solutions to the scalar wave equation may be viewed as two subproblems. 
First, the 0-matrix has to be created from a given dielectric structure, and secondly 
the eigenvalues, /3^, and the corresponding eigenvectors, have to be found. We 
have chosen to make our implementation in Matlab, because it provides effective tools 
for solving both these problems. 

To give a more precise description of our implementation, we consider an example 
where an arbitrary dielectric structure has been discretisized in a 100 x 100 grid. 
In this case becomes a matrix with 10000x10000 elements, which indicates that 
the finite difference approach is very demanding in terms of memory consumption. 
However, as most of the elements of are zero, it is advantageous to store as 
a sparse matrix, which is easily done with the Sparse-type in Matlab. For an n x 
n dielectric structure, we need to store elements in the full representation, while 
only ~ elements are required in the sparse representation. Obviously, this gives 
rise to a dramatic decrease in the memory consumption, and thereby a corresponding 
increase in the size of the dielectric structures that may be examined. 

The second problem we are faced with, concerns the search for eigenvalues, /3^, 
and their corresponding eigenvectors In the sparse representation, a complete 
diagonalization of the 0-matrix may be performed, but this straightforward method 
has the disadvantage of being numerically heavy (and thus time consuming) since 
it calculates all n'^ eigenvalues. Typically, we are only interested in solving for a 
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few eigenvalues, e.g the largest values of and this is a common feature of several 
advanced eigensolver libraries. In our implementation we use the EIGS command in 
Matlab which is based on the ARPACK-library [TJ. As a further advantage, the EIGS 
command has the option of finding eigenvalues around a specified value, which may 
be particularly useful once a region with guided modes has been found. 

3.2. Index-guiding Fibres 

As a first example, the finite-difference method is applied to an index-guiding 
micro-structured fibre. Figure \^ shows the square dielectric structure used in the 
calculation, which we have chosen to discretisize in a 128 x 128 grid. The dielectric 
structure has a width of 3\/3A, where A is the hole spacing, and it consists of air 
holes with a diameter of 0.4A placed in a silica background (n=1.45). A single air 
hole has been omitted in the center of the structure to create a waveguide core. In 
the case of index-guiding fibres, the fundamental mode corresponds to the eigenvector 
with the largest eigenvalue. Figure |2t shows the field distribution of the fundamental 
mode calculated at a normalized wavelength of A/A = 0.15, and it is seen to be well 
confined to the core region. 




Figure 2. (A) A dielectric structure of an index-guiding photonic crystal fibre, 
with a normalized holediameter of D/K = 0.4. For the calculation the structure 
has been discretized in 128 X 128 points. (B) The field distribution of the 
fundamental mode, calculated at a normalized wavelength of A/A = 0.15. 

In figure|31we have mapped out the effective mode index of the fundamental mode 
for several values of the normalized wavelength. For comparison, we have included 
a finite-difference calculation, where the width of the calculation domain has been 
increased to G-y/SA. Finally, we have included fuUy-vectorial results for a fully-periodic 
hole-structure, with a repeated defect, obtained in a plane-wave basis using periodic 
boundary conditions The dielectric structure used in this calculation consist of 6 
simple cells, and is thus similar to the structure in figure 

By comparing the finite-difference calculations for the small and the large 
calculation domain, we are able to see the effect of the truncation. For the 
small calculation domain, we find that the calculated value of the effective mode- 
index decreases rapidly for long normalized wavelengths. This is because the field 
distribution extends to the edge of the calculation domain, and thus the assumption of 
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Figure 3. Comparison of the mode-index for the fundamental mode in an index- 
guiding PCF with a hole diameter of 0.4A. The dotted curves are calculated 
by the finite-difference method (FDM) for two different widths of the dielectric 
structure, while the solid curve is calculated by the plane- wave method (PWM). 

a zero amplitude on the edge is no longer valid. Consequently, the field will penetrate 

into the air holes, and thereby cause a lowering of the effective mode index. By 
increasing the width of the calculation domain, we are able to shift the onset of this 
effect towards larger values of the normalized wavelength. 

In the comparison between our scalar finite-difference approach and a fuU- 
vectorial method, we find that the scalar approach gives qualitative correct results 
and accounts well for the overall wavelength dependence. However, especially for 
the long wavelengths the simple approach becomes inaccurate because of the scalar 
approximation, and the scalar approach is seen to overestimate the value of the 
effective mode-index. For A ^ A the strong proximity of the air-holes require a correct 
treatment of the air-silica boundary conditions which can only be quantitatively 
accounted for by a fuUy-vectorial approach. 

3.3. Photonic Bandgap- guiding Fibres 

In the literature, it is often argued that accurate modelling of photonic bandgap 
materials requires the use of a fuU-vectorial method. This is true for all photonic 
crystals of practical interest, because it is required that they have a large index- 
contrast in order to obtain large bandgaps. However, this may lead to the 
incorrect conclusion that photonic bandgaps and defect modes are pure full- vectorial 
phenomena. Prom a theoretical point of view bandgaps do not arise as a consequence 
of coupling between field components at a dielectric interface, as it is the case in a 
fuU-vectorial method. Rather, they are a fundamental property of applying the wave 
equation to a periodic waveguide structure, and thus photonic bandgaps and defect 
modes can also exist in a scalar calculation. 
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The question is how weU a scalar calculation actually approximates the results 
in photonic bandgap fibres, where the scalar wave equation is obviously not a correct 
representation of the actual physical problem. In order to examine this, we have 
chosen to apply our finite-difference method to the extreme case of airguiding photonic 
bandgap fibres. The dielectric structure used in our calculation is shown in figure 

It has a width of G-s/SA and consists of an air-silica cladding structure with a 
holediameter of I? = 0.88A. The core defect is made by inserting a large air hole with 
a diameter of 2.65A. For the calculation, we discretized the structure in 256 x 256 grid 
points. We have chosen this specific structure, because it is known to support guided 
modes in a full- vectorial calculation |12| . 




Figure 4. (A) The dielectric structure of an airguiding photonic bandgap fibre 
discretized in 256 X 256 points. The cladding structure has a hole diameter of 
0.88A while the core defect is created by inserting a large air hole with a diameter 
of 2.65A. (B) The field distribution of the fundamental mode in the airguiding 
photonic bandgap fibre, calculated at a normalized wavelength of A/A = 0.7. 

A disadvantage of this finite-difference implementation is that there is no simple 
way of calculating the position of the photonic bandgaps. Therefore, we have used 
a full-vectorial planewave-method to calculate the bandaps of an infinite triangular 
lattice with a hole diameter of 0.88A. Once the position of the photonic bandgaps 
have been located, it is possible to search for a defect mode. As already mentioned, 
the EIGS-comma.nd in Matlab has the useful ability to find eigenvectors around a 
specified eigenvalue. This is particularly useful for bandgap fibres, since the defect 
mode appears as an isolated eigenmode inside the boundaries of the photonic bandgap. 

By choosing a normalized wavelength of A/ A = 0.7 and searching for an eigenvalue 
for which (3 ^ k, the scalar finite-difference method finds a defect mode localized to 
the core region. The field distribution of this defect mode is shown in figure ^p. 
For simplicity we have chosen to depict the absolute value of the field distribution. 
Therefore, the 6 lobes surrounding the central defect do in fact have the opposite 
amplitude of the field inside the defect. These 6 resonances are a common feature of 
airguiding fibres, and they are also found in a full-vectorial calculation. 

In figure we have mapped out the effective mode-index for a range of the 
normalized wavelength. For comparison we have included the guided modes and the 
bandgap edges from a full-vectorial calculation. Both methods are found to predict 
the existence of a fundamental and a second-order mode, and a reasonable agreement 
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Figure 5. Comparison of the scalar finite-diflerence method (FD-method) and 
the full-vectorial planewave-method (PW-method) for an airguiding photonic 
bandgap fibre. A strong resemblance between the methods is found for both the 
fundamental and the second-order mode, but the finite-difl^erence is seen predict 
reasonably wider bandgaps. 

is found between the results of the two methods. However, we generally find that the 
finite-difference method overestimates the values of the effective mode-index. 

The major difFc'rc;iic;c; between the full-vectorial and the sc;alar calculation is that 
the latter is seen to predict significantly increased bandgaps. This naturally gives 
rise to a much wider range in which the structure supports a confined mode. The 
bandedges in the scalar calculation have been found as the modes that lie just above 
and just below the defect modes. As the bandedge-modes are infact cladding modes, 
and thus well distributed over the entire cross section of the dielectric structure, they 
are more affected by the truncation of the 0-matrix than the well confined defect 
modes. Wc have tried to increase the width of the calculation domain and also the 
number of grid points, but this is not found to have any influence on the overall result. 
It is therefore believed that the increased bandgap size is a consequence of the scalar 
approximation. 

The flnal result indicates, that although a scalar approach provides great insight 
to the basic physics of photonic bandgap fibres, it cannot reveal the complete picture. 
This is not really surprising. However, we still believe that this method is of great 
interest, mainly due to it simple implementation. Also, the model can easily be 
expanded to include periodic boundary conditions and a full-vectorial implementation 
is feasible as well. 

4. Conclusion 

The fleld of photonic crystal fibres has by now existed for almost a decade and the 
results of the research efforts will probably soon move inside the classroom and also 
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be found in text-book material on fibre optics and electro-magnetic theory of photonic 
crystals. This also calls for simple approaches to modelling of micro-structured optical 
fibres which are easy to implement and which without too much effort can produce 
results which are in qualitative agreement with those observed in real micro-structured 
fibres. We believe that the present work provides such a simple approach. 

In order to limit the complexity of the mathematical formulation of the problem, 
we have considered a scalar wave-equation which is solved by means of the most simple 
numerical approach to differential equations; the finite-difference approximation. With 
appropriate boundary conditions this results in a matrix eigenvalue problem. The 
matrix, 0, of this problem is highly sparse and has very simple analytical matrix 
elements which only depend on the lattice spacing, the frequency, and the dielectric 
function. Thus, no implementation of complicated basis functions is required. For 
a given frequency w the propagation constant /3 can easily be found by finding the 
eigenvalues of 0. By the aid of standard numerical routines, we are able to solve 
for a specific number of eigenvalues closest to a specified value. This is particularly 
useful for bandgap guiding fibres, where the modes of interest corresponds to either 
the smallest or the largest eigenvalue, but are placed in an interval determined by the 
photonic bandgap edges. 

In conclusion we find that the presented finite-difference method, apart from 
being simple to implement and despite the extremely simple model, is able to provide 
qualitative correct results for both index-guiding and photonic bandgap guiding fibres. 
The latter case is quite surprising, since modelling of photonic bandgap effects is 
normally associated with full- vectorial methods. Furthermore, we find that the method 
is robust and it is relatively simple to incorporate periodic boundary conditions, or to 
expand the model to a fuU-vectorial method. This holds interesting prospects for a 
future development of this method. 
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